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Abstract. We study the effects of energy transport in the Sun by asymmetric dark matter 
with momentum and velocity-dependent interactions, with an eye to solving the decade- 
old Solar Abundance Problem. We study effective theories where the dark matter-nucleon 
scattering cross-section goes as and with n = —1,0,1 or 2, where Urei is the dark 
matter-nucleon relative velocity and q is the momentum exchanged in the collision. Such 
cross-sections can arise generically as leading terms from the most basic nonstandard DM- 
quark operators. We employ a high-precision solar simulation code to study the impact on 
solar neutrino rates, the sound speed profile, convective zone depth, surface helium abun¬ 
dance and small frequency separations. We find that the majority of models that improve 
agreement with the observed sound speed profile and depth of the convection zone also reduce 
neutrino fluxes beyond the level that can be reasonably accommodated by measurement and 
theory errors. However, a few specific points in parameter space yield a significant overall 
improvement. A 3-5 GeV DM particle with a si oc (f' is particularly appealing, yielding 
more than a 6(T improvement with respect to standard solar models, while being allowed by 
direct detection and collider limits. We provide full analytical capture expressions for q- and 
Urei-dependent scattering, as well as complete likelihood tables for all models. 
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1 Introduction 

A non-relativistic relic dark matter particle, with a mass of a few GeV or more, is the 
leading candidate to explain astrophysical and cosmological phenomena ranging from cluster 
kinematics and galactic rotation curves, to gravitational lensing and the heights and positions 
of the acoustic peaks in the cosmic microwave background. Dark matter (DM) may have 
been produced in a similar way to Standard Model (SM) particles, either via chemical freeze- 
out (as in the weakly-interacting massive particle ~ WIMP - scenario) or via an initial 
asymmetry, analogous to baryogenesis (as in the asymmetric DM - ADM - scenario). If 
either of these scenarios is correct, it is possible that DM interacts weakly with SM particles. 
Such interactions would be seen most easily as a small elastic scattering cross-section between 
DM and quarks. 

The search for DM-quark interactions has been the focus of terrestrial direct detection 
experiments such as DAMA [1], CoGeNT [2], CRESST-II [3] and CDMS II [4], who have all 
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reported excess events above their expected backgrounds. On the other hand, XENONIO [5], 
XENONIOO [6], COUPP [7], SIMPLE [8], LUX [9] and SuperCDMS [10] have all established 
strong limits on the DM-nucleon cross-section, seemingly in contradiction to the excesses 
observed by other experiments. These underground detectors are typically most sensitive to 
DM particles with masses of order = 50 to 100 GeV, which lead to the largest recoil 
energies on the heavy nuclei used as targets. For the same reason, they are also best suited 
to probing fast-moving DM particles, resulting in a threshold velocity of a few tens ofkms”^ 
for an incoming DM particle to create a nuclear recoil event. 

Collisions between DM and nuclei can also lead to capture and accumulation of DM 
in the solar core. For this to occur, collisions between DM and nuclei in the Sun must 
result in sufficient energy transfer for the DM velocity to be brought below the local escape 
velocity. The kinematic region probed by the Sun is quite different to the one probed by direct 
detection: optimal energy transfer, leading to optimal capture rates, occurs for DM particle 
masses closely matching the solar composition, i.e. a few GeV. Because DM gains speed as 
it falls into the solar potential well, the low-velocity tail of the local DM distribution is the 
dominant contributor to solar capture; the opposite is true for direct detection experiments. 
Direct detection and solar physics are therefore highly complementary laboratories for the 
study of DM-quark scattering. 

1.1 Dark matter in the Sun 

The effect of DM on stars has been the subject of investigation for some time (see [11-13] 
for reviews). There are two main scenarios where DM might create observable effects on the 
Sun: 1) the capture and annihilation of WIMP-like particles; and 2) the accumulation of 
an “asymmetric” species. The first case is characterised mainly by a search for high-energy 
neutrinos with Ej, ~ [14-31] and MeV-scale decay products [32, 33], using detectors such 

as IceCube and SuperKamiokande. Annihilation of DM in a star can also provide an energy 
source through the release of other SM particles [11, 34-47], leading to changes in the core 
structure. The capture rate of DM in the Sun is however far too low for the energy released 
this way to have any significant effect on solar structure [11]. 

Our focus in this paper is therefore on the case where DM either cannot annihilate, 
or does so far more slowly than it is captured. In this case DM is effectively asymmetric 
(i.e. ADM, [48, 49]), and large quantities of DM may have built up inside the Sun over its 
lifetime. The weak interactions of DM with quarks give the DM particles relatively large 
inter-scattering distances, making them potentially significant conductors of energy. Just 
like the capture process, energy transfer is most efficient for lighter DM masses (1 — 10GeV), 
as momentum transfer is maximised when the masses of the colliding particles are equal, and 
the Sun is mostly H (A = 1) and He (A = 4). It is interesting to note that this is roughly 
the mass range expected in the most common models of ADM [50], in order to explain the 
1:5 cosmological relic baryon-to-DM density ratio. In contrast, earth-based direct detection 
experiments lose sensitivity in this range, as they make use of high-mass elements such as 
germanium (A ~ 72) and xenon (A ~ 131), for which high incoming DM velocities are 
necessary to create a measurable recoil. 

There is a window of elastic DM-nucleon scattering cross-sections [51-55] at u ~ 10“^® 
cm^ for spin-independent (SI) couplings (10“^^cm^ in the spin-dependent - SD - case) for 
which a is large enough to allow sizeable capture of DM, but small enough that the mean 
interscattering distance is still large. A large inter-scattering distance allows heat to be 
redistributed away from the solar core by DM-nucleon scattering. This has the effect of 
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reducing the temperature of the solar core Tc, and increasing the central density and pressure. 
Changes of state variables near the solar centre affect the production of neutrinos from 
fusion processes, especially the flux of neutrinos, which goes as Tc with f5 ~ 20-25. The 
temperature and density changes in the core reduce the local sound speed, and force an 
overall mass redistribution that impacts the solar structure at other radii. This leads to 
modifications of the sound speed profile over the entire depth profile of the Sun, and shifts 
the height of the base of the convection zone. Both the sound speed profile and the depth of 
the convection zone have been independently measured using helioseismology. 

It has been shown that high-precision solar evolution models including capture, trans¬ 
port and (minimal) annihilation of DM can be built to satisfy the observed solar age, radius, 
and luminosity [53-55]. At the same time, it appears possible for the inclusion of DM in 
such models to affect the (less well constrained) ®B flux in an observable way, and even 
improve agreement with the observed sound speed profile and the depth of the convection 
zone. Neither of these latter two observables are well reproduced by standard solar models 
computed with the latest surface compositions [56-61]. This issue, known as the “solar com¬ 
position” or “solar abundance” problem, has been brought about by the 20-25% reduction 
in the measured solar metallicity in recent years [62-73], and is one of our motivations for 
this paper. 

Although several studies have indicated that ADM can alleviate some of this tension, 
the cross-sections required to do so are typically far higher than allowed by limits from direct 
detection. Here we investigate whether broader consideration of the kinematic structure of 
the DM-nucleon vertex might provide a way around this. In the process, we provide first 
rigorous limits from solar physics on such DM models, which we refer to as ‘generalised 
form factor dark matter’. In a separate paper [74] we discussed a specific realisation of 
momentum-dependent dark matter that leads to a Qa improvement over the Standard Solar 
Model (SSM). We revisit this model in Sec. 6.1. 

1.2 Generalised form factor dark matter 

The kinematic differences between direct detection and the Sun become even more marked 
if the DM-nucleon interaction is not assumed to be independent of the DM-nucleon relative 
velocity, or the momentum transferred in the collision. There is indeed no guarantee that 
the standard SI and SD operators correctly represent the DM-quark interaction. In particle 
physics, the interaction cross-section generally depends on the centre of mass energy and the 
transferred momentum, parameterised using the Lorentz-invariant Mandelstam variables s, t 
and u. In the non-relativistic limit, these become the centre of mass momentum, proportional 
to the relative velocity Ureh and the transferred momentum q = Ap. As these are small 
quantities, the constant term usually dominates in a series expansion of the cross-section. 
However, many models with non-trivial dependencies on Urei and q exist, typically motivated 
by theoretical arguments or attempts to reconcile experimental results. 

To make quantitative predictions, a specihc form of a{vre\, q) must be chosen. Because 
we wish to remain as general as possible, we choose to focus on couplings of the form a oc 
and a oc q^"-, with n = { — 1,1, 2}. The and forms are respectively called p-wave and 
d-wave interactions, and correspond to the cases where the initial state particles possess 1 and 
2 units of relative angular momentum, respectively. These are always present, but normally 
only dominate when all lower-order terms in the scattering matrix element - including the 
constant (s-wave) term - are suppressed due to cancellations. Cross-sections depending on q 
can arise, for example, from a non-zero particle radius (the analogue of a nuclear form factor), 
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from parity-violating couplings like XlbXQQi Xlii'lbXQx^Q and X^Q■, or from a 
small anapole or dipole interaction between the dark and visible sectors [75-87]. We refer to 
the class of models where DM-nucleon scattering cross-sections depend on some combination 
of q and/or Uj-ei as ‘generalised form factor DM’ because it generalises the effects of form 
factors to arbitrary powers of q and Urei- 
Concretely, we focus on the couplings 


a = ao 



(1.1a) 


and 



(1.1b) 


These can lead to either spin-dependent (SD) or spin-independent (SI) interactions, depend¬ 
ing on the axial structure of the DM-nucleon interaction vertex. The normalisation ctq must 
be defined with respect to some reference velocity vq or momentum go- We will choose 
vq = 220kms“^, the typical halo DM velocity, and go = 40 MeV, corresponding to a nuclear 
recoil energy of around 10 keV in an underground direct detection experiment. 

The DM.-nucleus cross-section is related to the above DM-nucleon cross-sections via: 


0'N,i 


m2(m^ + "inuc)^ 


O'SI^i + o-SD 


4( + 1) 

SJi 




( 1 . 2 ) 


where Ai and Jj are respectively the atomic number and total angular momentum of nuclear 
species i\ (Sp^i) and {Sn,i) are the spin expectation values of its proton and neutron systems. 
Given the dependence of Eq. 1.2 we will find that, in spite of the Sun’s small metallic- 
ity, spin-independent DM can have a significantly larger effect than a spin-dependent DM 
candidate, which couples mostly to hydrogen.^ Beyond a few studies, e.g. [88], this fact has 
not been emphasised very much in the literature. The full effect of a momentum-dependent 
cross-section will furthermore depend crucially on the composition of heavier elements. 

The full impact of momentum and velocity-dependent DM on solar observables has not 
been studied before. The authors of Ref. [30] computed the effect of such couplings on the cap¬ 
ture rate of spin-dependent DM, and computed neutrino fluxes from an annihilating species. 
However, they did not include a treatment of heavier elements, nor of the crucial energy trans¬ 
port by DM once captured. In Refs. [89, 90], the authors computed capture and transport 
rates for ADM models with long-range interactions. To account for the impacts of the non¬ 
trivial scaling of these cross-sections with momentum and velocity, they employed effective 
cross-section scaling factors to decouple the cross-sections entering the capture and transport 
calculations, and avoid modifying the standard velocity-and-momentum-independent treat¬ 
ment of capture and energy transport. As we show later, it happens that this rescaling can 
indeed be done without any loss of generality for capture in the velocity-dependent case, but 
it is not possible in the momentum-dependent case. It is also not possible to account for the 
effects on energy transport of either a velocity or momentum dependence in this manner; 
rather, a full recalculation of the transport coefficients must be performed [91]. 

The structure of this paper is as follows: in Section 2, we review the capture equations 
for DM in the Sun, and present the necessary modifications for velocity and momentum- 
dependent scattering of DM with nucleons. In Section 3 we review the theory of conductive 
heat transport by DM developed in Refs. [91, 92] , along with its application to solar modelling. 
Section 4 describes the DarkStec computer code that we have developed for simulating the 
effects of generalised form factor DM on the Sun. We present results in Section 5, and discuss 

^See also caveats to this treatment in Sec. 2.2.1. 
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their implications with regards to the Solar Abundance Problem, current experimental limits 
and previous work on the topic in Section 6. We summarise in Section 7. 


2 Capture of dark matter by the Sun 


2.1 Standard (velocity and momentum independent) treatment 

The population of DM particles in the Sun is follows the differential equation 

= C^(t) - A(t) - E(t), (2.1) 


where is the capture rate, A(t) is rate at which annihilations occur and E{t) represents 

evaporation. Unless DM is strongly self-interacting (not the case we consider here, but 
discussed in Ref. [93]), (^©(t) does not depend on the DM population already captured 
by the Sun. A(t) is the rate of annihilation, and is proportional to the square of the DM 
population. Here we consider the case where DM is fully asymmetric, so A{t) = 0, although 
we comment briefly in Section 4.1 on the implications of allowing DM to self-annihilate. 

Evaporation occurs when a DM particle gains enough energy from a scattering event to 
overcome the Sun’s gravitational potential and escape, so E{t) is linear in the DM population 
(being simply the product of the single-particle evaporation probability and the number of 
candidates for evaporation). Evaporation requires a significant gain in momentum, as the 
typical velocity of a thermalised DM particle is ~100kms“^, whereas the escape velocity in 
the solar core approaches 1400 km s“^. This means that evaporation is only significant if DM 
is similar in mass to the nucleus with which it scatters in the Sun. In practice this means that 
DM about the mass of helium (4 GeV) and lighter is typically most prone to evaporation, but 
other masses closely matched to significant elements in the Sun (C, N, O, Ee ~ 12, 

14, 16, 56 GeV) can also in principle be affected [94]. The evaporation rate depends on the 
interaction cross-section, mean free path and thermal regime (LTE vs Knudsen); extending 
the standard analyses to generalised form factor DM is therefore non-trivial. In practice 
the rate of evaporation is extremely low in almost all cases where the nuclear scattering 
cross-section is allowed by direct detection; for very specific analyses it should be taken into 
account, but for the purposes of this paper we assume that the evaporation rate is zero. We 
intend to return to this point in detail in a future paper. 

We now turn to the capture rate as it is implemented in our simulations. As we are 
using the DarkStars code [95], we closely follow Refs. [11, 96]. We take the local distribution 
function f{u) of DM to be Maxwell-Boltzmann, with dispersion uq = 270kms“^. In the 
frame of the Sun, moving at uq = 220kms“^ relative to the Galactic rest frame, 


/ 3 \ 3/2 / 3(uq -|- u^)\ sinh(3nnQ/uQ) 

\2y TT^/^m^UQ 2uq J 3uuq/uq 


( 2 . 2 ) 


As it falls into the gravitational potential well of the Sun, a DM particle acquires a velocity 
w = It becomes gravitationally captured by the Sun if it loses enough 

kinetic energy in a scattering event for w to fall below the local escape velocity Vescii^, t). Eor 
this to occur, the fractional energy lost by the DM particle A = 2E-{i/, corresponding 
to nuclear recoil energy Hr, must be in the interval 


w'^ ~ m^w'^ ~ p?j^ ’ 


(2.3) 
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where /x = m^^/mnuc and /x± = (/x ± l)/2. The local capture rate of particles with velocity 
w is the sum over nuclear species at radius r, of the rate of scattering in the interval of Eq. 
2.3, i.e. 


Q{w) 


2 

m^w 


E 


t) 




' m-^u^ 12 


Fi(ER)|2dER, 


(2.4) 


where |Ti(ER)p is the nuclear form factor. For hydrogen this is a constant, whereas for 
heavier nuclei we use the usual Helm form factor 


Fi(ER)|2 = exp 


(2.5) 


Here Ei is a constant quantity for each nuclear species x, given by 

( 2 . 

Integrating over the phase space, the total capture rate of DM in the Sun is then 

CQ{t) = A-k [ wi}{w) du dr. (2.7) 

Jo Jo u 

For a constant cross-section and a Maxwell-Boltzmann velocity distribution, this can be 
solved analytically [11, 96]. However, in the following we generalise the above equations to 
allow momentum and velocity-dependent usi and crsD-^ In this case, numerical integration 
becomes necessary. 


Ei = 


5.8407 X 10 


-2 


mAr,j(0.91m]^j -h 0.3)^ 


GeV. 


2.2 Velocity and momentum-dependent treatment 

We begin with a momentum-dependent cross-section of the form Eq. 1.1a, with a oc To 
include such a cross-section in Eq. 2.4, the constant cross-section must be replaced with Eq. 
1.1a and the dependence on the nuclear recoil energy moved inside the form-factor integral. 
This explicitly illustrates the equivalence of a momentum-dependent cross-section to a change 
in form factor, and we refer to the corresponding integral and associated multiplicative factors 
as the ‘generalised form factor integral’ (GFFI). For hydrogen, |E(Er)P = 1, so the change 
required is 


GFFIn=o = I \FiEn)\^dEn 


G'EE/„^o,h 


VgoJ 2/x" 



dA, 


( 2 . 8 ) 


where we have expressed the form factor integral in terms of A instead of Er for compactness. 
We do this by noting that the transferred momentum is q = so the fractional 

energy change can be written A = fiq'^/p'^, where p = m^w is DM particle’s incoming 
momentum. Performing the integral in Eq. 2.8 yields 



^We will use cr as a shorthand for either ctsi or ctsd, as the spin and kinematic dependence can be factorised. 
For an explicit treatment, see [97]. 
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Figure 1. Left: effective enhancement/suppression of DM capture in a standard model of the Sun 
(AGSS09ph; [61]), due to a momentum-dependent cross-section, as computed with Eq. 2.11. The 
black line at = 1 corresponds to tr = const. Solid lines are spin-independent (SI) and dashed lines 
are spin-dependent (SD), coupling only to hydrogen. Right: effective enhancement/suppression of 
energy transport as defined in Eq. 3.10 , in the LTE regime, with cto = 10“^® cm^. We note that 
away from LTE, the behaviours reverse, and a q~^ cross-section actually yields an enhancement with 
respect to the constant case. This can be seen in Fig. 2. The full effect of a g-dependent cross-section 
is then the combined effect of the left and right-hand panels. 


For heavier elements, the integrand includes the Helm form factor (Eq. 2.5), so 


GFFIn^0,ijLH 


pV" Ej 
qo) {B^iY 


r 


(^ + n,B 



-r 


^1 -b n, H 



( 2 . 10 ) 


where B = ^rriyw'^/Ei, and r(m,x) is the (upper) incomplete gamma function. To gain a 
more intuitive understanding of Eqs. 2.9 and 2.10 we show in Eig. 1 the overall enhancement 
or suppression to the capture rate with respect to the constant case. 




capture 


Z,nA^GFFIn=oy 


( 2 . 11 ) 


for three values of n yf: 0. In this example we just take w = w{r = 0) and use the present-day 
AGSSOOph Standard Solar Model^ [61, 70]. Here fi represents the fractional composition of 
each species with atomic number Hj. Eor the spin-dependent cross-sections, the sum over 
the species i only includes hydrogen. Of course the actual capture rate must be accurately 
integrated over the entire star, and we compute it precisely with Eq. 2.7 in our simulation. 
In each case, the overall behaviour roughly follows the (p/(?o)^"' dependence, as heavier DM 
will have gained more momentum as it falls into the solar gravitational well. 

In the velocity-dependent case (Eq. 1.1b), where a oc the modihcation to the 
capture rate is much simpler. The partial capture rate Q,{w) oc cr is simply transformed by 
the replacement 

( 2 . 12 ) 


^Publicly available at http://www.ice.csic.es/personal/aldos/Solar_Models.html. 
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The integral Eq. 2.7 can then be evaluated to obtain the modified capture rate. We later 
show with our solar simulation code that the enhancement due to velocity-dependence agrees 
very well with what is obtained using the standard (n = 0 ) treatment and an average of the 
cross-section throughout the volume of the star 

cr{vrei) ^ W{wo)) = / a{wo)pQ{r)r‘^ dr, (2.13) 

^"0 Jo 

where wo = is the typical velocity acquired by an in-falling WIMP as it reaches 

radius r. This indicates that a cross-section proportional to or - typically thought of 
as a suppression - actually yields a considerable enhancement. We note that cr(urei) has been 
erroneously approximated to a{{wo)) rather than {a{wo)) in the literature [90], although 
[uesc(?' = 0 )/uo]^" does at least provide the correct enhancement to the capture rate to within 
a factor of a few. 


2.2.1 Form Factors 

We finally comment on the use of the Helm form factor 2.5 in our capture equations. Very re¬ 
cent work [97] has shown that a more accurate computation of the nuclear response functions 
to DM-nucleus scattering can have a substantial impact on the capture rate. These response 
functions amount to corrections to the Helm form factor by terms proportional to powers of 
y = Eii/2Ei. These will dominate at large momentum transfers. We have checked that in 
most cases considered here - and in all cases where the effect of DM is large enough to modify 
the solar observables - this quantity is much smaller than one. This is because the low DM 
masses that yield the largest conduction effects are mainly sensitive to the leading, constant 
term in the NR expression. However, we caution the reader that accurate modification of 
the capture rate for DM with larger masses than those considered here, should make use of 
the full nuclear response functions detailed in Ref. [97]. 

We have also used the standard approach for spin-dependent capture in the Sun, wherein 
it is assumed that hydrogen is the dominant contributor to the capture rate. It was shown in 
[97] that, while this is strictly true for velocity-dependent cross sections (non-relativistic oper¬ 
ator O 7 ) and holds to a few percent accuracy for a constant cross section (O 4 ), momentum- 
dependent interactions can lead to a different behaviour. In the case (Oq), nitrogen 
dominates capture by an order of magnitude or more for all DM masses. The net effect is 
that the spin-dependent capture rates computed here for have been underestimated with 
respect to the case when the full nuclear response functions are included. 


2.3 The geometric limit 


In all cases, the total effective cross-section of the Sun to collisions with DM particles cannot 
exceed the geometric “cutoff” 

o-max = 7rR|(f), (2.14) 

which corresponds to the case where the Sun is optically thick to DM. This places a funda¬ 
mental limit on the capture rate 


C*max(f') — TtR, 


l-c 

Ut) 

Jo 


feiu) 


u 


w‘^{u, Rq) du 


(2.15) 


3 rriy. 


' _3 “© 

e 


6 OCnVI©-|-R0(uq-|-3uq) 

-uo H-z--^Erf 

vr RqUq 


3^ 

2 uo 
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The actual capture rate to be used must therefore be the lesser of Eqs. 2.15 and 2.7. Assuming 
a steady radius and a local DM density of py- ~ 0.38 GeV cm“^, the maximum amount of 
DM that can accumulate in the Sun by its current age is therefore 


iVx,max 1.5 X 10^^ 



~ 1.3 X 10"^^ 



where A^b — Mq/ttIp is the total number of baryons in the Sun. 


(2.16) 


3 Conductive energy transport by dark matter 


If enough DM is captured by the Sun, its large typical inter-scattering distance means that 
it is a more efficient carrier of heat over long distances than ordinary baryonic material, so it 
can act as an additional mechanism for heat transport alongside photons. The formalism we 
use to compute the effect of microscopic energy transport by DM conduction was developed 
by Gould and Raffelt [92]. By solving a perturbative expansion of the Boltzmann collision 
equation (BCE) in the Sun’s gravitational potential, they showed that the thermal conduction 
by a weakly-inter acting species can be expressed in terms of two quantities: a dimensionless 
molecular diffusivity a{p), and thermal conductivity k{p), where p = m^/mnuc is again the 
ratio between the DM and nucleon masses. If more than one nuclear species is present, a 
and K get replaced with effective values, which are weighted by the number densities of each 
species in the plasma mixture at each height in the Sun. 

In the local thermal equilibrium (LTE) regime, where Ixi^) is much smaller than both 
the inverse of the local temperature gradient |VlnT(r)| and the DM scale height r^, the 
values of a and k are found by solving the first order expansion of the BCE. This is done 
in terms of the quantity e = /^(r)|Vlnr(r)|, via the formal inversion of the Boltzmann 
collision operator C{u,r), where C{u,r)F{u,r) represents the net change in the DM phase 
space distribution F{u, r) due to collisions with nuclei. The collision operator is defined 
via phase space integrals over collision rates, meaning that the dependence of the collisional 
cross-section cr on Urei and q must be explicitly included in the calculation. In Ref. [92] Gould 
and Raffelt computed and tabulated a{p) and k{p) for a constant scattering cross-section 
(n = 0). In Ref. [91], we extended the formalism to include velocity and momentum- 
dependent cross-sections (n 7 ^ 0 ), showing that the resulting changes in a, k and can have 
potentially large effects on heat transfer in the Sun. 

The equilibrium distribution of DM particles in the gravitational potential (/>(r) of the 
sun is given by [11, 91, 92] 


^X,LTe(^) — ''^x,LTe(0) 


rT(r )1 

3/2 

[t(o)J 

exp 




dr' 


d0(r^) 


kBTir') 


where r = 0 represents the centre of the Sun. The conductive luminosity is: 
Ax,lte(?’) = 47rr'^C^'^{r)K{r)nx,LTE{r)lx{r) 


kBT(r)Y^^ , dT(r) 
Kb- 




dr 


(3.1) 


(3.2) 


where the factor accounts for a velocity-dependent [C = ro/uT(r)] or momentum-dependent 
[C = Qo/'n^x'’^T{r)] cross-section, and r'r(r) = yj2kBT{r)/rrix is related to the typical thermal 
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Figure 2. Illustration of the transition from the local thermal equilibrium (LTE) to the Knudsen 
(non-local, isothermal) regime of energy transport by DM scattering, for momentum-dependent (left) 
and velocity-dependent scattering {right). Leftwards of the peak in each curve corresponds to the 
Knudsen regime, whereas rightwards is the LTE regime. The total energy transport is plotted for a 
fixed DM mass {m^ = lOGeV) and number ratio of DM to baryons {n^jn^, = 10“^^). Solid lines are 
spin-independent couplings, whereas dashed lines represent the spin-dependent case where the DM 
scatters only on hydrogen. This has the effect of increasing the mean free path, leading to a transition 
to the Knudsen regime at a much higher value of (Tq- 


velocity [92]; vq and go are respectively the reference velocity and momentum defined in Eq. 
1.1. The rate of energy transported per unit mass of stellar material is: 


ex,LTE (J’) 


1 

4:Trr^p{r) 


dT^,LTE(?^) 

dr 


(3.3) 


This quantity is usually expressed in units of ergs g“^ s“^. 

If the condition l^. <C Vy. is violated, LTE is no longer valid and the system exists in the 
Knudsen regime of non-local transport, where the DM distribution is essentially isothermal. 
It was shown in Refs. [92, 98] by Monte Carlo simulation that Eqs. 3.1 and 3.2 can be corrected 
to properly account for a large Knudsen number K = Ix/^x- Here we have formally defined 


Following [11, 52], we dehne: 


/ 3fcBr, 

\27rGpcmxJ 


W) 

m 



(3.4) 


(3.5) 

(3.6) 


where Kq = 0.4 and r = 0.5 are empirical values taken from the numerical results of [98]. 
The DM distribution is then a combination of the isothermal and LTE distributions: 


%(t) = f(-?^)rax,LTE + [1 - f(-fS^)] %,iso, 


(3.7) 
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where 


n. 


' 7 ^ 


lx,iso{r,t) - 

Finally, the Knudsen-corrected luminosity is: 


(3.8) 


total (r,t) = f(iF)f)(r,t)L;^,LTE(r,t). 


(3.9) 


This should then be used in place of T^^lte in Eq- 3.3 to compute the energy injected 
or removed at each radius by DM-nucleon collisions. In Fig. 1 we illustrate the effective 
enhancement or suppression of transport in a toy solar model due to a momentum-dependent 
cross-section, plotting 


_ f \e{r,n ^ 0)\r‘^dr 

-'transport = “Wj - ——j , 

J |e(r, n = U)|r^ar 


(3.10) 


for the three cases of momentum-dependent cross-section, with do = 10“^® cm^. Once again, 
this is illustrated using a present-day SSM with the AGSS09ph abundances [61, 70]. Note 
that this cross-section leads mainly to transport near the LTE regime. As do is decreased 
further towards the Knudsen regime, the behaviour reverses, as the enhancement provided by 
the longer inter-scattering distance is overcome by the Knudsen suppression. We illustrate 
the Knudsen behaviour in Fig. 2 by plotting the total energy transport for different types of 
generalised form factor DM as a function of the cross-section do, for a constant DM-to-baryon 
ratio n^/ub = 10“^^. This shows how the peak energy transport varies in each model. We 
note that this peak occurs for a much smaller cross-section in the spin-independent case, due 
to the reduced mean free path caused by scattering with helium and metals rather than just 
hydrogen. 

The full effect of ADM on energy transport is then a combination of the effects illustrated 
in the left and right panels of Fig. 1, keeping in mind the degree of non-locality indicated by 
Fig. 2. 


4 The DarkStec solar dark matter code 

In order to accurately model the effects of generalised form factor DM on solar observables, 
we implemented full velocity- and momentum-dependent DM capture and energy transport 
in the high-precision solar evolution code GARSTEC [99, 100]. We took DM routines from 
the public dark stellar evolution code DarkStars [95], producing a hybrid code DarkStec. 

GARSTEC is the descendant of the legendary Kippenhahn code. Numerical aspects and 
physics inputs are described in detail in Ref. [99] and the modified version of GARSTEC used 
for this work is the same described in Ref. [100]. Here we just give a summary of the most 
relevant physical inputs. It includes the nuclear energy generation routine exportenergy.f^, 
updated with the astrophysical factors recommended in the Solar Fusion II [101] compilation. 
It makes use of the Opacity Project radiative opacities [102], complemented at low temper¬ 
atures with those from Ref. [103]. The equation of state is the 2005 update of OPAL [104]. 
Microscopic diffusion of elements, including gravitational settling, thermal and concentration 
mixing, is treated according to Ref. [105]. 

The calibration of a solar model generally implies adjusting a number of free parameters 
in the model to match an equal number of observables. In the present case, the observables 

^Publicly available at http: //www. sns. ias. edu/~ jnb. 
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are the present-day (r© = 4.57 Gyr) solar luminosity Lq, solar radius Rq and the metal- 
to-hydrogen mass fraction {Z/X)q. The latter is a critical quantity, as it determines the 
composition, namely the metallicity, of the calibrated model. In this work, we adopt the 
photospheric solar abundances from Ref. [70], for which {Z/X)q = 0.0180. The free param¬ 
eters in the model are the mixing length parameter omlt and initial helium and metal mass 
fractions lini and Zi^i respectively. The latter two suffice to determine the initial abundances 
of all elements in the model because the relative metal abundances are taken from Ref. [70] 
with Zini acting as the normalisation factor, and Xini -|- lini + -^ini = 1 by definition (Xini 
being the initial hydrogen abundance). 

In practice, the solar model calibration starts with a homogeneously-mixed 1 Mq pre- 
main sequence model that is evolved assuming no mass loss until it reaches tq, the current 
solar system age. At that age, model predictions are compared with the observables, and a 
Newton-Raphson scheme is implemented to iteratively find the solution. This is generally 
achieved to 1 part in 10® within two to three iterations for standard solar models (i.e. with no 
DM). More iterations are necessary as the effects of DM become more important. In the most 
extreme cases, no physical solutions are found (e.g. resulting in negative Tmi)- Note that each 
iteration requires four evolutionary calculations to evaluate the partial derivatives needed for 
the Newton-Raphson scheme. Each evolutionary calculation requires 700-800 timesteps. At 
each timestep, the solar structure is discretized in about 2000 shells. These requirements for 
the integration of solar models, both in spatial and time resolution, guarantee a numerical 
precision better than 1% in all model predictions [58]. 

DarkStars [95] is a Fortran95 package that implements capture, annihilation and energy 
transport by regular (n = 0) WIMP dark matter in a general stellar evolution code, as 
described in Refs. [II, 39, 40, 106]. The capture routines were originally adapted from 
DarkSUSY [107]. The underlying evolutionary code [108] is a FortranQO rewrite of the the 
venerable Fortran77 Cambridge STARS package [109-111]. These codes use the relaxation 
method to solve the coupled ID ordinary differential equations of stellar structure over an 
adaptive grid. DarkStars is the state of the art in dark stellar evolution for the n = 0 case, 
as it features the full capture calculation (Eq. 2.7) for SI and SD scattering on the 22 most 
important elements, various options for the DM velocity distribution (including user-defined 
distributions), and proper Gould-Raffelt treatment of conductive energy transport (Eq. 3.2), 
including self-consistent density profiles and the Knudsen-dependent interpolation between 
the LTE and isothermal (non-local) regimes (Eqs. 3.7, 3.9). 

For DarkStec, we adapted the DarkStars capture routines to implement capture of gen¬ 
eralised form factor DM (Eq. 2.9-2.10, 2.12) rather than just the n = 0 case. We used the 
conductive transport routines from DarkStars essentially unaltered, except that we included 
the additional factor of in Eq. 3.2 and utilised the a and n tables that we computed 
earlier for n ^ 0 [91]. 

At each regular GARSTEC timestep, DarkStec computes the total DM capture rate by 
solving Eq. 2.7, assuming a local halo DM density of 0.38 GeV/cm®. This input rate is then 
used to update the DM population in the Sun following Eq. 2.1. DarkStec uses the numerical 
version of the capture routines from DarkStars to evaluate the modified capture equation Eq. 
2.7. It computes the thermal diffusivity and conductivity coefficients a{r,t) and n{r,t) at 
each height in the star by interpolating in the tables of Ref. [91], as these quantities depend 
on the specific mixture of plasma species with which the DM particles interact. The code 
then computes the DM density n^(r, t) using a{r,t) and K{r,t), which it uses to determine 
energy transport. DarkStec then interpolates the resulting values of (Eq. 3.3) to the 
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grid used by GARSTEC, where they are treated as an additional energy source at each height 
in the star. 

We considered seven ADM models with SI interactions with nucleons, and seven with 
SD interactions: the constant cross-section case, a oc and a oc with n = {—1,1, 2}. 
For each coupling, we simulated one solar model for each point on a grid of and a^. We 
computed a subset of models using a stringent convergence criterion of one part change in 
10^ for the the solar luminosity, radius and surface metallicity {Z/X), with 4kyr and lOMyr 
minimum and maximum time steps. Because it was extremely computationally expensive to 
do every simulation including the full treatment of capture and transport at this accuracy, 
we carried out all other simulations with a more relaxed convergence criterion of a part in 
10^, with minimum and maximum time steps of 40kyr and 40Myr. We then corrected these 
lower-accuracy results to consistently reproduce the observables and chi-squared values of the 
higher-accuracy models, using the systematic differences we saw between models computed 
both ways. In total, our calculations took ~ 1.5 CPU years. 

4.1 Annihilation 

DM models with momentum- and velocity-dependent nuclear scattering need not necessarily 
be entirely asymmetric. To investigate the implications of energy injection from annihilation 
in such models, we also implemented annihilation in DarkStec, following [11] and DarkStars. 
In general however, allowing annihilation simply weakens the limits that we obtain from solar 
physics on generalised form factor DM, and does not improve the overall fit to solar data in 
any significant way. We therefore do not show these results, nor discuss annihilation beyond 
this subsection. It is worth noting that although we assume zero annihilation cross-section 
for all the results we show, some small (sub-thermal) annihilation cross-section is certainly 
still allowed in each model, and would make no impact on the solar observables. 

5 Limits on generalised form factor dark matter from solar physics 

In this section we present a systematic overview of the results obtained from our simulations. 
In each subsection, we review the effects of velocity and momentum-dependent asymmetric 
dark matter on a specific observable. These are: the boron-8 and beryllium-7 neutrino fluxes, 
the depth of the convection zone rcz, the sound speed profile Cs(r), the small frequency 
separations and the surface helium abundance Is. Agreement between the predicted and 
observed values of the neutrino fluxes is typically unchanged or worsened by thermal transport 
by DM. The same is true for Is. In contrast, the predictions of rczj the sound speed profile 
and the small frequency separations are often closer to the observed values when conductive 
energy transport by DM is included. In Sec. 5.7 we construct a combined likelihood, which 
encompasses the overall improvement or degradation in the fit to solar data for each cross- 
section form and combination of uo and m^. We also present a few interesting benchmark 
cases in which the agreement is significantly improved. 

For every form of the cross-section, we ran simulations over a grid of DM masses and 
cross-sections: = 5, 10, 15, 20 and 25 GeV, and each decade in uo between 10“^° and 

10“^*^ cm^. The colour scales in the figures of this section are interpolations from this grid. 
For some specific cases, where low-mass points gave a good overall improvement over the 
Standard Solar Model, we extended the mass axis down to 1 GeV. We once again caution that 
although evaporation is expected to have an effect near these small masses, a full kinematic 
analysis would be necessary to determine its exact and uo dependence for each model; 
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Figure 3. Cross-section normalisation ao at which dark matter capture saturates the geometric 
limit of Eq. 2.15 for momentum-dependent (left) and velocity-dependent (right) cross-sections. Spin- 
independent couplings are shown with solid lines, and spin-dependent with dashed lines. The SI q~^, 
and cases are not shown, as they saturate the capture rate below ag = 10“^° cm^ (the smallest 
cross-section that we simulated). We note that by using the standard assumption that spin-dependent 
DM is only captured by hydrogen, we have underestimated the capture rate for g^; the suppression 
to capture could therefore be as much as an order of magnitude smaller in this specific case. 


the importance of evaporation for these models remains essentially unexplored. Given the 
importance of kinematic matching with individual nuclei, we caution against placing too 
much store in quick estimates of this effect. 

We show the impacts of generalised form factor DM on capture rates and their saturation 
in Fig. 3, and on observables in Figs. 4 to 17. 

Simulations that led to a reduction in <p'^ of approximately ~60% or more did not 
converge, because the effects on the overall solar evolution are too large for the numerical 
solver to handle. In many cases, the excess energy transport actually led to a density inversion 
in the core. In some of these models, a solution probably exists, and could be found with a 
better solver. In others, the Sun probably cannot exist as a stable body. We simply mask 
the the non-converged region in blue in all our plots. 

5.1 Saturation of capture 

A velocity or momentum dependence in the nuclear scattering cross-section has a striking 
impact on the rate of DM capture by the Sun. The minimal ag required to render the Sun 
opaque to dark matter - and thus to saturate the capture rate Eq. 2.15 - can be lowered by 
as much as two orders of magnitude for a cross-section, and approximately one order of 
magnitude in the q~‘^ and cases. On the other hand, negative powers of Urei and positive 
powers of a spin-dependent g-coupling yield a large suppression in the capture rate, seen as a 
much larger required cross-section to achieve saturation. We illustrate the required value of 
iTo in Fig. 3, based on the output of our simulations. Interestingly, a spin-independent or 
q'^ cross-section can yield an enhancement or suppression of the capture rate depending on 
the DM mass; this is a consequence of the behaviour illustrated in Fig. 1, and simply reflects 
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the fact that g'-dependent scattering is most efficient for large momentum transfers, which is 
much easier when the DM mass is closely matched with the masses of heavier elements. 

5.2 Solar neutrino fluxes 

In general, energy transport by DM removes energy from the solar core and reduces its 
temperature, causing a reduction in neutrino production rates. Given its strong temperature- 
dependence, the neutrino production rate is the first place to look for changes in solar 
observables. 

In Figs. 4 and 5, we show the effect on the ®B neutrino flux for spin-independent and 
spin-dependent dark matter with velocity and momentum-dependent couplings. Fluxes are 
normalized to the observed value, (j^'^ohs ~ ^ cm“^s“^; lighter colouring represents 

a reduced flux. Although the measurement error on the ®B neutrino flux is only 3% [112], 
the overall uncertainty from modelling is around 14%. We add the absolute uncertainties in 
quadrature to obtain the la region. White lines therefore represent the la isocontour where 
the neutrino flux falls below ~ 85% of the measured value; black lines are the 2a (71%) 
contours. 

Although the flux of ^Be neutrinos is not as temperature-sensitive as that of ^B neutri¬ 
nos, they can also be used as a weaker, independent, probe of the solar core temperature. In 
Figs. 6 and 7 we show the ratio of the predicted ^Be neutrino fluxes to the observed value 
4>Beohs ~ ^ cm“^s“^. Here again we plot la and 2a contours using the theoretical 

and observational uncertainties, which are respectively 7% and 5% for ^Be neutrinos. 
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Figure 4. The ratio of the predicted ®B neutrino flux to the measured value obs = ^ x 10® 
for each type of spin-independent dark matter coupling defined in Eq. 1.1. In every case the white and 
black lines show the isocontours where the flux is respectively 1 and 2cr lower than the observed values, 
based on observational (3%) and modelling (14%) errors, added in quadrature. The cross-sections 
are normalized such that a = or cr = cFQ{qlqQY'^ , with uq = 220kms“^ and go = 40MeV. 

Simulations carried out in the masked regions did not converge, due to the significant heat conduction 
by the DM particles, leading in extreme cases to density inversions in the core. 
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Figure 5. As per Fig. 4, but for spin-dependent couplings. We note that by neglecting capture 
on elements other than hydrogen in the SD case, we have underestimated the capture rate (and net 
effect) for the case. See discussion in Sec. 2.2.1 
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Figure 6. The ratio of the predicted ^Be neutrino flux to the measured value = 4.82 x 

10® for each type of spin-Independent dark matter coupling defined in Eq. 1.1. In every case 

the white and black lines show the isocontours where the flux is respectively 1 and 2a lower than the 
observed values, based on observational (5%) and modelling (7%) errors, added in quadrature. The 
cross-sections are normalized such that a = ao{v/vo)^'^ or cr = o'o{q/qo)^^, with vq = 220kms“^ and 
qo = 40 MeV. 
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Figure 7. As per Fig. 6, but for spin-dependent couplings. 
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As expected, the effect of ADM on neutrino fluxes is larger for low DM masses, mainly 
due to the increased efficiency of momentum transfer between DM and H/He as the DM takes 
on a similar mass to these two dominant nuclei. Ignoring evaporation, transport becomes 
even more efficient as mass is decreased to 1 GeV for SD dark matter; in the SI case, the 
largest effect is seen at ~ 3 to 4 GeV, around the average nucleus mass in the Sun. 

The dependence on the cross-section is also as expected: below a minimum value of 
uo, not enough DM can be captured to significantly affect the production rate. For 
constant cross-sections, this is around 10“^® cm^ in the SI case, and 10“^® cm^ for SD dark 
matter. This highlights the importance of heavier elements, to which only the SI DM couples 
in our simulations: the large fraction of helium, combined with the dependence of the 
DM-nucleus cross-section (Eq. 1.2) means that the behaviour of DM inside a star is highly 
sensitive to elements heavier than hydrogen. 

The upper edges of the regions shown in Figs. 4-7 highlight the transport behaviour: if 
the cross-section falls below a certain value, it allows the DM to more efficiently penetrate 
the dense plasma, carrying energy from hot to cooler regions. The location of the maximum 
reflects a combination of being near the transport efficiency peak (where the transition from 
the LTE to the Knudsen regime occurs) while simultaneously maintaining a large capture 
rate through a large enough cross-section. 

The impacts of velocity and momentum-dependent scattering of dark matter can be 
understood in terms of this balance between capture and transport. Positive powers of Urei 
lead to an enhancement in the capture rate, moving the window in which DM transport 
has a significant effect on neutrino fluxes to lower cross-sections. However, the important 
suppression in the overall transport rate, as seen in Fig. 5 of [91], yields an overall suppression 
in the effect relative to the constant case. The opposite behaviour is seen for a oc v~J, 
although the enhanced transport is not sufficient to compensate for the ~10“^ suppression 
in capture. 

Again, the effect of a momentum-dependent cross-section is more subtle. Although 
the boosted capture rate for a q~^ cross-section moves the region of interest to lower cross- 
sections, closer to what is allowed by underground experiments, the overall suppression in 
transport means that there is very little overall effect on solar observables. Positive powers 
of q fare much better, with an enhancement both in capture and transport. 

Comparison of Figs. 4 and 5 highlights the importance of heavier elements in such 
processes. In most cases, they lead to a significant enhancement of the capture and transport 
rate; however, we note in comparing the q~^ plots that they can also inhibit energy transport, 
by reducing the mean free path for conduction. 

5.3 Depth of the convection zone 

The boundary between the convective and radiative zones rcz is the location at which the ra¬ 
diative and adiabatic temperature gradients are equal. Above this location, the temperature 
gradient becomes super-adiabatic, hydrostatic equilibrium breaks down and convection sets 
in. The amount of energy transported by DM is negligible at heights r ~ rcz- However, the 
changes in the density and temperature gradient near the core lead to a small but significant 
increase in the radiative temperature gradient at much larger radii, causing it to exceed the 
adiabatic gradient at a slightly lower depth and shift the lower boundary of the convection 
zone downwards. 

In Figs. 8 and 9, we show the ratio of the location of the lower boundary of the convection 
zone predicted by each model to the value inferred from helioseismology, rcz,© = 0.713 i?©. 
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Figure 8. Ratio between the modelled and measured location of the bottom of the convection zone 
rcz^ for spin-independent couplings. Darker regions represent a better ht to the observed value than 
the Standard Solar Model. The white and black lines represent the contours at which the predicted 
value falls within Itr and 2a of the measured value, respectively. The theoretical uncertainty on rcz 
(0.004i? 0 ) is much larger than the experimental error ( 0.001 i?©), so the former dominates when we 
add them in quadrature. 
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As per Fig. 8, but for spin-dependent couplings. 
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Although the error on the helioseismological inference is just 0.001 Rq, the theoretical error on 
the modelled rcz,© is much larger: crcz,0,th = 0.004i?©. Adding these errors in quadrature, 
we plot in white and black the contours containing the regions where rcz,© is within Icr and 
2 (7, respectively, of the inferred value. 

As can be seen in Figs. 8 and 9, the difference between the depth of the convection zone 
in the Standard Solar Model (rcz,© = 0.722 Rq) and the depth inferred from helioseismology 
amounts to almost a 3(7 discrepancy. Except for a small region at high cross-section 
(do = 10“^^ cm^), spin-dependent DM struggles to bring rcz to within much better than 2a 
of the inferred value. However, energy conduction from SI interactions does substantially 
better: and models can produce good agreement with the observed value of rcz,©) 

leading in some regions to agreement at better than Icj. 

5.4 Surface helium abundance 

The surface helium abundance Yg is another observable that has fallen into disagreement with 
observations since the revision of the solar composition. Our SSM prediction is = 0.2356, 
whereas the observed value is 0.2485 ±0.0034. With theoretical errors of ±0.0035 taken 
into account, this amounts to a discrepancy of 3a. The addition of dark matter does very 
little to change Yg, and for most models the discrepancy remains at the 3a level. For some 
very few cases (not shown), it actually worsens the discrepancy by up to an additional ~2(7. 
This only occurs for models where the fit to sound speed observables is also substantially 
worsened, notably for large SD and q^ cross-sections. The reduction in Yg is caused by a 
corresponding increase in Xi, the initial hydrogen fraction. The increase in Xi is demanded 
by the reduction of the core temperature, which requires a greater amount of hydrogen to be 
present in the core in order to maintain the same nuclear burning rate as in the SSM, and 
to thereby match the observed solar luminosity Lq. 

Because it changes little, we do not show the contour plots for the surface He abundance. 
For completeness though, we include it anyway in our full likelihood computation in Sec. 5.7. 


5.5 Sound speed profile 

In Fig. 10 we show the deviation of the modelled radial sound speed of some models that 
produced the best overall fits with respect to the measured values from helioseismology. To 
illustrate the agreement of each model with the observed sound speed profile, we construct 
an effective chi-squared measure 


Xc„ 


E 


ri 


[Cs,model(^2) Cg^hel. (^i)] 


cr. 


^s,hel. 


(d) 


(5.1) 


where the sum goes over 5 equally-spaced radial points between r = 0.1 Rq and 0.67 Rq. We 
do not include radii below 0.1 i?© because the reconstructed sound speed in this region is 
highly uncertain due to the low number of low-degree (low-/) p-modes reaching the innermost 
radii of the solar core [113]. The upper limit of the range is the location of the largest 
discrepancy in the sound speed of the Standard Solar Model, at the base of the convection 
zone; above this point essentially all models agree very well with the observed sound speed 
because the temperature gradient is adiabatic and therefore does not depend on the detailed 
composition of the Sun. We added errors from modelling and inversions for each point in 
quadrature. We obtained modelling errors by using models for which one input parameter 
was varied at a time to obtain partial derivatives at each radial point, and then combined 
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Figure 10. Deviation of radial sound speed profile from values inferred from helioseismology. Left: 
the best overall fits with constant spin-dependent (SD) and spin-independent (SI) cross-sections. 
Right: best-fit models of the three couplings returning the best overall p-values (see Tab. 1): {q^, 
SD): = 5 GeV, cto = 10“^° cm^; (g^, SI): = 3 GeV, ctq = lO”^"^ cm^; and (u^, SI): = 5 

GeV, CTO = cm2. 


quadratically given the uncertainty in the AGSS09 abundances and errors in each parameter 
quoted in [114]. The errors on inversions were taken from [115]. The values of this meant 
to show the relative improvement to the sound speed modelling, are shown in Figs. 11 and 
12 . 

From these figures, along with the sound speed profiles illustrated in Fig. 10, it is clear 
that the addition of momentum or velocity-dependent dark matter to the solar model does 
indeed help alleviate the discrepancy between modelling and observation in some specific 
cases. We will return to these cases in Sec. 6.1. 


5.6 Frequency separation ratios 

The inverted sound speed profile Cs{r) obtained from helioseismological measurements is not 
very accurate near the solar core because not enough low-Z p-modes are available for a precise 
and accurate inversion. Instead, information about the solar core can be gained by using the 
so-called frequency separation ratios. A large advantage of using these ratios is that, unlike 
individual frequencies, they are not affected by the detailed structure of the solar surface, 
which is poorly described by solar models [116]. This is because for radial order n 3> 1 and 
low angular degree I, the surface effects are functions of the eigenfrequency, so they cancel out 
when considering frequency differences between modes of similar frequencies. In addition, 
by taking ratios of appropriate frequency differences, the solar core structure becomes the 
dominant effect in the observed signal [117, 118]. In particular, two very useful quantities 
are the frequency separation ratios 


where 
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Figure 11. Effective sound-speed Xc defined in Eq. 5.1, showing the improvement in goodness-of-fit 
of the radial sound speed profile in solar models with different types of spin-independent ADM. The 
value of the best fit point is indicated by a green line on the colour bar and is shown as a green 
star in each figure. Ax^ = 2.3 and 6.18 contours, corresponding to Icr and 2cr deviations from the 
best fit, are respectively shown in white and black. 
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Figure 12. Same as Fig. 11, but for spin-dependent couplings. 
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Figure 13. Small frequency separations ro 2 (top) and ri 3 (bottom) as defined in Eq. 5.2, for the best 
fit constant (left) and generalised form-factor dark matter models (right). The latter correspond to 
the best-fit models of the three couplings returning the best overall p-values (see Tab. 1). Predictions 
are compared with helioseismological observations from the BiSON experiment [117]. Inner black 
error bars correspond to observational error, whereas outer (green) bars also include modelling error. 
Below each figure we show the residuals with respect to BiSON data, in units of the total error. 
The rriy^ = 3 GeV, case with ctq = 10“^^ cm^ yields the best improvement, bringing the largest 
discrepancy from nearly dcr to little more than a standard deviation. 
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and A/(n) = The sound speed gradient is weighted by 1/r in the integral above, 

so ro 2 (n) and ri 3 (n) are most sensitive to changes there. 

We show a set of these ratios in Fig. 13 for the same examples shown in Fig. 10. These 
are compared with the values of ro 2 (n.) and ri^^n) measured by the BiSON experiment [118]. 
Model uncertainties are computed following the same procedure as for the sound speed. As 
can be seen in the residual subplots, the SSM overestimates each of these ratios by as much as 
4(7, once modelling errors are included. Thermal transport by ADM can smooth the sonnd 
speed gradient, yielding smaller differences in Eq. 5.3. As in the case of the sound speed 
profile, the best improvement is provided by scattering with rriy^ = 3 GeV, do = 10“^^ 
cm^. In this model, the largest discrepancy in ro 2 and ri 3 falls to barely la. 

In Figs. 14 and 15, we show the overall ability of different models to fit the observed 
frequency separations. We quantify this with the combined chi-squared Xro 2 T Tria) where 




2 

re,e +2 


E 


[ri,i+2,th.(n) - r£,£+2,obs.(n)]^ 


(5.4) 


Formally, there is a correlation between values of ro 2 (u) and r 13 ( 71 ) with common eigenfre- 
quencies. In practice, however, the correlation is minimal (the Pearson correlation coefficient 
is always lower than 0.1), so we assnme they are independent when computing Xro 2 + Aris- 
Figs. 14 and 15 show that all SI and SD couplings do exhibit areas of parameter space where 
the small frequency separations can be brought into better agreement with observations than 
in the SSM. In many cases however, these are not the same parameter combinations preferred 
by other observables. The notable exception to this is the SI coupling, which also clearly 
produces the best agreement with the observed frequency separations. 


5.7 Combined limits 

It is clear from the results we have presented in Secs. 5.2-5.6 that the disagreement between 
model predictions and helioseismology can be ameliorated, at the cost of a slightly worse 
agreement with neutrino fluxes. Surface helium is also affected to a lesser extent. It is 
therefore instructive to construct an overall likelihood function based on the fit of these data 
by our models. We define the combined chi-squared as: 

2 (1 “'/'B,obs/‘i^B)^ , (1 “‘?^Be,obs/'i^Be)^ 

T = -2-+-2- 

^B ^Be 

(rcz - ■rcz.obs)^ (Ts - Ys^obs)^ 

^CZ ^Ys 

+Xro2+Tri3- (5-5) 

where ~ x 10®cm“^s“^, 4>Beobs ~ ^ 10®cm“^s“^, rcz.obs = 0.713i?© and 

Ys,obs = 0.2485. The uncertainties include both the observational and modelling errors, 
added in quadrature. These are given individually in the last two lines of Table 1. We 
include the small frequency separations r£^i^ 2 , but not the sound speed profile, as the latter 
is less precise and is correlated with the former. 

We show the combined limits from Eq. 5.5 in Eigs. 16 and 17, comparing with limits 
from direct detection where such results exist (i.e. in the SI case; [30]). We discuss these 
results, and the comparison with other limits, in more detail in the following section. 


- 28 - 








cTsi cx const. 



Figure 14. Combined likelihood Xro 2 ^^13 the small frequency separations dehned in Eq. 5.2, 
for different models of ADM with spin-independent scattering cross sections. The green line on the 
colour bar shows the value of the best ht, indicated by a green star in each panel, whereas white, 
black and cyan lines show the preferred regions at 1 , 2 and 4 a conhdence levels, respectively. At 
low masses, many DM models can yield improvements in the overall chi-squared with respect to the 
Standard Solar Model, but the model provides a clear best Ht. 


- 29 - 









































(To (cm^) (To (cm^) cq (cm 


(TsD oc const. 




Figure 15. Same as Fig. 14, but for spin-dependent couplings. 
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Figure 16. Full combined likelihood (Eq. 5.5) incorporating ®B and ^Be neutrino flux measurements, 
the surface helium abundance Fsi the radius of the convection zone and the small frequency separa¬ 
tions. The value of the best fit point is shown as a green star, and indicated by a green line on the 
colour bar. = 2.3, 6.18 and 19.33 contours corresponding to Icr, 2cr and 4cr deviations from the 
best fit are shown in white, black and cyan, respectively. 


- 31 - 








































(To (cm^) (To (cm^) cq (cm 





Figure 17. Same as Fig. 16, but for spin-dependent couplings. 
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Table 1. Standard Solar Model (SSM) and best fit (b.f.) values for each of the models we consider, 
along with observable quantities. The DM mass and cross-sections are in GeV and cm^, respectively. 
The boron-8 neutrino flux c/ig is in units of 10“® cm“^ s“^ and the beryllium-7 neutrino flux 
is expressed in 10“® cm“^ s“^. The full chi-squared is defined in Eq. 5.5 and includes the neutrino 
fluxes, surface helium abundance Ig, depth of the convection zone and small frequency separations. 


Model 

{nix, *^ 0 ) 6 ./. 

</>B 

^Be 

Rcz/Rq 

W 


P 

SSM 

- 

4.95 

4.71 

0.722 

0.2356 

287.8 

< 10-^'J 


const. 

(15,10-^0 

3.48 

4.37 

0.721 

0.2348 

122.2 

< 10-^'J 



(3,10-38) 

4.51 

4.59 

0.721 

0.2350 

144.1 

< 10-1° 



(3, 10-37) 

3.78 

4.29 

0.718 

0.2327 

27.5 

0.85 

O'SI 


(10, 10-31) 

3.29 

4.27 

0.721 

0.2348 

98.2 

1.1 X lO-i" 


v~‘^ 

(5,10-34) 

3.91 

4.44 

0.721 

0.2347 

86.4 

5.0 X 10-° 



(5, 10-38) 

3.31 

4.22 

0.719 

0.2337 

135.0 

< 10-1° 



(5,10-38) 

4.30 

4.54 

0.721 

0.2349 

114.9 

3.6 X 10-1° 


const. 

(5,10-33) 

3.36 

4.27 

0.720 

0.2341 

100.2 

5.8 X 10-8 



(5,10-38) 

4.43 

4.59 

0.722 

0.2353 

174.5 

< 10-1° 


q^ 

(5,10-30) 

3.64 

4.35 

0.720 

0.2343 

73.1 

2.5 X 10-^ 

crSD 

q^ 

(15,10-30) 

3.81 

4.47 

0.721 

0.2351 

140.7 

< 10-1° 



(10,10-32) 

3.39 

4.33 

0.721 

0.2347 

114.5 

4.2 X 10-1° 



(5,10-33) 

3.75 

4.39 

0.721 

0.2344 

72.4 

3.1 X 10-^ 



(5,10-37) 

4.80 

4.68 

0.722 

0.2355 

238.7 

< 10-1° 

Obs. 

- 

5.00 

4.82 

0.713 

0.2485 

- 

- 

Obs. error 

- 

3 % 

5% 

0.001 

0.0034 

- 

- 

Model error 

- 

14% 

7% 

0.004 

0.0035 

- 

- 


In Table 1, we give a break-down of the values of each observable at the best-fit parame¬ 
ters of each model, along with the best-ht and implied p-value. The best-fit regular SI and 
SD cases shown in Table 1 differ from the best-fit cases in Ref. [74] because here we actually 
choose them on the basis of the full likelihood Eq. 5.5. In Ref. [74] we only computed the 
small-frequency separations for models that provided the best fits to all other observables. 
Because the fits to different observables are inconsistent in the models with constant cross- 
sections, including the frequency separations shifts the best-fit masses and cross-sections. 
This is to be compared to the SI case, where adding in the small frequencies merely makes 
the best-fit point even better. 

We also provide as supplementary material online a complete table of the boron and 
beryllium neutrino fluxes, surface helium, convective zone radius and small frequency sepa¬ 
rations for all of our models, in addition to the partial and total chi-squared values defined 
in Eq. 5.5. This table can be used for quick lookup and interpolation for, e.g. global fits of 
DM models. 

6 Discussion 

6.1 Solving the Solar Abundance Problem 

Our results of Section 5 indicate that the impacts of different types of generalised form 
factor DM on solar observables is quite varied, and not always straightforward. What is 
clear, however, from Eigs. 16 and 17 and Table 1, is that some models do indeed lead to a 
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much better overall fit than the Standard Solar Model alone. Although there is always some 
reduction in the goodness of fit to observed neutrino fluxes, this can be accommodated in 
many cases by the theoretical errors. 

In Fig. 10 we showed the reduction of the discrepancy between simulated and measured 
sound speed profiles of the best-fit models of the couplings that give the best overall improve¬ 
ment compared to the Standard Solar Model. On the left, we showed the best models for 
constant spin-dependent and spin-independent cross-sections, and on the right we showed the 
best-fit models of the three best generalised form factor couplings, according to the likelihood 
defined in Sec. 5.7. We find that the best improvement comes from a scattering cross-section 
that is proportional to the square of the transferred momentum q. This is confirmed to 
very high significance by the small frequency separations, seen in Fig. 14. At low masses 
~ 3 —5 GeV, an asymmetric DM candidate with a = (To((?/40MeV)^ and (Tq = 10“^^ cm^ 
leads to a substantial improvement with respect to the SSM. The best fit occurs at 3 GeV, 
and leads to a convective zone depth rcz = O.TlSi?©, only 1.2a away from the inferred 
value from helioseismology. The and ^Be fluxes are respectively 3.78 x 10® cm“^s“^ and 
4.29 X 10® cm“^s“^, both within 2a of the measured values. Surface helium remains low, at 
Yg = 0.233. The large improvement to the fit to Cs(r) can be seen in the right-hand panel 
of Fig. 10, and to the small frequency separations probing the solar core in Fig. 13. This 
specific case yields an overall 6 a improvement compared to the Standard Solar Model, as we 
discussed explicitly in Ref. [74] . 

In the next section, we will discuss constraints from other sectors, and show that the 
best-fit model does indeed constitute a viable candidate for reconciling solar modelling, 
helioseismology and spectroscopy. It is worth noting, however, that the masked regions of 
non-convergence in our contour plots may hide parameter combinations that yield even better 
fits than this; careful work on the numerical solver is required to explore this possibility. 

Other g-dependent models, as well as constant and velocity-dependent models we exam¬ 
ined, nearly all led to poor overall fits: the effect on and Yg is always too large, in spite of 
the improved fit to helioseismological measurements. Every other ADM model we examined 
led to some significant improvement over the SSM, but failed to provide p-values greater than 
10“^. As pointed out in Ref. [90], a v~J SD cross-section does indeed lead to an improved 
fit to the sound speed profile and rcz for a small region in mass and cross-section space (see 
Figs. 9 and 12). Unfortunately, this case is disfavoured by our overall likelihood, mainly 
due to excessive reduction in the core temperature, which is reflected in a greater-than 3a 
reduction in the ®B neutrino fluxes, as well as significant reduction in ^Be neutrinos. In fact, 
we find that a SD cross-section yields a much more plausible overall fit to observables 
than SD scattering, resulting in a p-value of 3.1 xlO ^ as compared to 4.2 xlO (vs. 
0.85 for the SI model). It is also important to note that the region where SD improves 
helioseismological agreement is very close to the boundary beyond which our models failed 
to converge. This issue, also seen by the authors of Ref. [54], could indeed be hiding a better 
overall fit that is simply beyond the reach of our solver. Given the great reduction in neutrino 
fluxes that we observe in the bordering area however, we remain somewhat sceptical of this 
possibility. 

6.2 Comparison with collider and direct limits 

The dashed curves in Fig. 16 show independent limits from direct detection on SI scattering, 
as derived in Ref. [30]. To the best of our knowledge, similarly model-independent limits on 
q and Urei-dependent SD couplings from direct detection are not available in the literature. 
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For momentum-dependent couplings, the direct limits suggest that all interesting masses 
(from the point of view of solar physics) above ~5 GeV are in tension with direct searches. 
Curiously though, the areas that lead to the best overall improvements in solar observables {q^ 
scattering, ~ 3-5 GeV, usi ~ cm^) are some of the few that survive the direct 

search limits. For velocity-dependent and constant couplings, the curves lie almost entirely 
below the plotted regions, indicating that in these cases, the entire section of parameter space 
that turns out to have an impact on solar physics is disfavoured. 

It is worth remembering, however, that limits from direct detection and solar physics 
probe very different parts of the dark matter halo velocity distribution: solar capture prefer¬ 
entially occurs from the low-velocity end of the distribution, whereas direct detection probes 
the high-velocity tail. The compatibility of solar and direct searches is therefore highly sen¬ 
sitive to the chosen velocity distribution, and might be modified substantially by choosing 
something other than the standard Maxwell-Boltzmann distribution. 

Limits on both momentum-dependent and velocity-dependent DM-quark couplings also 
exist from colliders [119-122], although in this case a specific effective operator or UV- 
complete model must be assumed in order to calculate definite event rates. In the clas¬ 
sification scheme of [123]([119]), operators F2(D2) and FIO(DIO) have leading contributions 
to the nuclear recoil rate from q'^ SI terms, operator F6(D6) from SI terms, operators 
F3(D3), F6(D6) and S2(C2/R2) from q^ SD terms, operator F4(D4) from q^ SD terms, and 
operator S4(G4) from SD terms. For the most part, ATLAS and CMS have not focused 
on such operators due to the weakness of the limits on them from direct detection. Limits 
from Tevatron data have been calculated by CDF [124] on the combination of F4(D4) and the 
‘regular’ constant SI scalar coupling XXQQ (Fl/Dl), but so far only a few phenomenological 
groups [119, 122] have placed limits on the other operators using up-to-date collider data. 

From the point of view of solar physics, given the preference we see for SI scattering 
in ameliorating the Solar Abundance Problem in Fig. 16, operators F2(D2) and FIO(DIO) 
would seem to be the most interesting. Existing limits on FIO(DIO) from ATLAS data [122] 
are quite stringent at niy. below a few hundred GeV, excluding mass scales of below a TeV 
for the physics associated with the particle mediating the interaction. Such limits would 
likely make achieving our preferred reference cross-section rather difficult. Limits on F2(D2) 
are much weaker though, with only mass scales below 42 GeV excluded at < 50 GeV 
(assuming that the contact operator approximation even remains valid in this regime). It 
therefore seems feasible that a DM particle interacting with quarks by this effective coupling 
(xI^XQQ) might explain the recent anomalies seen in some direct detection experiments, and 
help rectify the known discrepancies between theory and data in solar physics. 

6.3 Comparison with previous work 

The effect of DM heat conduction in the Sun, especially by asymmetric dark matter, has been 
studied in several works. Most recently, the main qualitative features of such models were 
pointed out by [125] while accurate computations using solar models were done in Refs. [53- 
55], with somewhat contradictory conclusions. In this section we will endeavour to settle the 
differences. Among these studies, only [54] used the proper thermal conduction mechanism 
of Gould and Raffelt [92]. Our energy transport and neutrino flux computations are 
in good agreement with that work. However, Ref [54] did not include the (small) effect 
of molecular settling, nor did they account for the uncertainty in the initial helium mass 
fraction Y, metallicity Z or mixing length parameter a, all of which are free parameters of 
solar models. Allowing Y, Z and a to vary has the effect of partially compensating for the 
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increased luminosity from DM transport, thus reducing the overall effect of DM and leading 
us to weaker overall constraints. 

In contrast, Refs. [53, 55] used the approach developed by Spergel and Press (SP, [126]). 
However, comparison with Monte Carlo simulations found this approach to be flawed [98]. 
This is due to the violation of all three assumptions that go into the SP approach: 1) a 
Maxwellian velocity distribution; 2) spatial uniformity of the WIMP velocity distribution; 
and 3) local isotropy. To further quantify the errors introduced by the use of the SP treatment, 
we have implemented the SP transport mechanism into our code, following Ref. [51] for a 
constant, spin-dependent (SD) cross-section. We plot the total energy transport for a given 
DM population in the Sun as a function of the spin-dependent cross-section in the left- 
hand panel of Fig. 18. Solid lines represent the accurate Gould and Raffelt (GR) approach 
(which we use in this paper), whereas dashed lines are the SP calculation. As in Fig. 2, the 
peak energy transport occurs at the transition between local (LTE) and non-local (Knudsen) 
transport. Although both approaches roughly yield the same location for this peak, SP 
slightly underestimates transport in the LTE regime, whereas it signihcantly overpredicts 
transport in the non-local regime. In the large cross-section LTE regime, the SP case predicts 
DM luminosities that are 1.5 (for my_ = 5 GeV) to 4 (m^ = 80GeV) times too small; in 
the Knudsen limit, this becomes an overestimate by a factor of 3 (m^ = 80 GeV) to 9 
{niy- = 5 GeV). It is worth noting that the largest discrepancies happen very near the location 
of interest, where low masses and cross-sections potentially lead to effects in the Sun and 
simultaneously evade constraints from direct detection. In the right-hand panel of Fig. 18 
we show the resulting discrepancy between predicted neutrino fluxes with SP and GR 
transport, as computed with our solar model. The red region is where SP overpredicts an 
effect, whereas blue corresponds to the LTE regime where SP underestimates the energy 
conduction. 

Two other recent papers [89, 90] have also investigated the impacts of specific dark 
matter models with non-constant cross-sections on solar physics. In both these papers, 
the authors employed existing capture and transport schemes (SP in [89] and GR in [90]) 
without properly accounting for the dependence of the capture and transport processes on the 
kinematic structure of the cross-section. Ref. [89] examined a magnetic dipole dark matter 
model with a cross-section that scales as relative to the standard constant case.^ 

They then applied the ‘momentum transfer cross-section’, which amounts to multiplying the 
cross-section by (1 — cos 9) in order to cancel the divergence induced by forward scattering 
when integrating over all scattering angles to obtain a total cross section. This process is 
only useful for obtaining a rough scattering rate, as it simply transmutes a momentum- 
dependence into an equivalent Urei-dependence (see e.g. calculations of k for q~‘^ couplings 
in [91] for details). The correct treatment is of course to include the full momentum- and 
velocity-dependent differential cross-section in the integral over the DM velocity distribution 
function. In the case of Ref. [89], this treatment simply converted a cross-section to a 

= Constant cross-section, effectively ignoring the momentum and velocity-dependence 
of the model. Similarly in Ref. [90], where the application of the momentum transfer cross- 
section treatment converted what was actually a cross-section in the long-range limit 

into a cross-section, relative to the standard case. The authors then took a 

thermal average over the velocity distribution and applied it as a constant scaling to the 
standard capture and transport treatments. As we have shown previously [91], this is valid 

®Here we remind readers that even the ‘constant’ case has Aa/dq^ oc i.e. whenever we refer to a q" or 
scaling, we mean relative to this familiar case. 
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Figure 18. Comparison between the Gould and RafFelt (GR) transport formalism and the incorrect 
Spergel and Press (SP) formalism, for a constant, spin-dependent cross-section. Left: Total energy 
transported by asymmetric dark matter from GR (solid lines) and SP (dashed lines) for a constant 
ratio of DM to baryons n^jn-b = 10“^® in a standard solar model. The SP formalism, which was 
shown to be flawed by [92, 98], overpredicts the transported energy by as much as a factor of 9 in the 
non-local (Knudsen) regime, while underestimating it in the LTE case. Right: comparison between 
predicted ®B fluxes once these transport mechanisms are implemented into the full DarkStec solar 
capture code. In the red region, using SP yields a larger reduction in the neutrino flux than GR. 


for the capture rate (or would be, if the cross-section were truly dependent only on Urei and 
not q), but qualitatively incorrect when dealing with conductive energy transport. 

In light of these shortcomings in the way capture and transport were modelled, and in 
the theoretical treatment of the scattering cross-sections involved, it is difficult to estimate 
how the results of Refs. [89, 90] would change were the kinematics of the cross-sections 
treated correctly. Calculating the transfer coefficients as per Ref. [91] and carrying out the 
full phenomenological analysis as we have here, for and couplings, would 

therefore be an interesting topic for future study. 

7 Conclusions 

In this paper, we have explored the prospects for solving the Solar Abundance Problem using 
generalised form factor dark matter, confirming that it can indeed provide a solution to this 
long-standing issue [74]. In the process, we derived the full capture rates in the Sun for dark 
matter with momentum- or velocity-dependent interactions with nucleons, and used solar 
physics to place limits on those couplings. This is the first extensive, physically-consistent 
analysis of the effects on the Sun of momentum- or velocity- dependent dark matter, taking 
into account the full capture kinematics, along with the self-consistent thermal conduction 
framework developed in Ref. [91]. We incorporated this framework into a new state-of-the 
art dark solar evolution code, DarkStec. 

Our best-fit model, asymmetric DM of ~ 3 GeV and ao = 10“^^ cm^ interacting 
with nucleons proportionally to q'^, shows excellent promise for solving the Solar Abundance 
Problem. Our model provides a remarkable ht, and a full investigation of the effects of 
evaporation, which has not been done before for non-constant cross-sections, should serve 
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as a cross-check on the consistency of such a scenario. The prospects for detecting such a 
particle at the next round of direct detection experiments, and in the next run of the LHC, 
also appear very promising. 
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